#include "header.h"

static int 		n, i, j, k, myid, np, ikx, iky, ikz;
static ptrdiff_t 	l_0_start;
static double		l, dx, height, width, m, kx, ky, kz, c, kk;
static int		flag, gFlag;
static MPI_Status	status;
static double		offset;

void init_ic(ctrl_ptr ctrl)
{
  n  = ctrl->n; gFlag = ctrl->gFlag;
  l  = ctrl->l; c = (2*PI/l);
  dx = l/n; m = ctrl->m;
  height = ctrl->b;
  width  = ctrl->a;
  offset = ctrl->offset;
  MPI_Barrier(COMM);  
}
void set_ic(ctrl_ptr ctrl, work_ptr in)
{

  init_ic(ctrl);
  double  *x = malloc((in->l_n0)*sizeof(double));
  double  *y = malloc(n*sizeof(double));
  double  *z = malloc(n*sizeof(double));
  double  rr;

  MPI_Comm_size(COMM, &np);
  MPI_Comm_rank(COMM, &myid);
  MPI_Barrier(COMM);
  //printf("Local n0 %d\n",in->l_n0);
  if (myid != 0) MPI_Recv(&flag, 1, INT, myid-1, myid, COMM, &status);
  for (i = 0; i < (in->l_n0); ++i )
  {
    x[i] = (myid*(in->l_n0)+i+0.5 + offset)*dx - l/2;
    //printf("myid:%d\t x[%d] = %f\tmyid*l_n0\n",myid,i,x[i]);
    for (j = 0; j < n; j++)
    {
      y[j] = (j+0.5+offset)*dx - l/2;
      for (k = 0; k < n; k++)
      {
         ikx = i + myid*(in->l_n0);
	 iky = j;
	 ikz = k;
 	 if (ikx > n/2) ikx = ikx - n;
	 if (iky > n/2) iky = iky - n;
	 if (ikz > n/2) ikz = ikz - n;
 	 kx = c*ikx; ky = c*iky; kz = c*ikz;
	 z[k] = (k+0.5+offset)*dx - l/2;
         rr = sqrt(x[i]*x[i]+y[j]*y[j]+z[k]*z[k]);
	 kk = kx*kx+ky*ky+kz*kz;
	 (in->wrk[k+n*(j+n*i)]) = height*cexp(-width*rr*rr);
	 if (gFlag == 0) (in->grv[k+n*(j+n*i)]) = cexp(-m*rr)/rr;
	 else (in->grv[k+n*(j+n*i)]) = 4*PI/(dx*dx*dx*(kk+m*m));
         //(in->grv2[k+n*(j+n*i)]) = 4*PI/(dx*dx*dx*(kk+m*m));
      }
    }
  }
  if ((myid == 0)&&(gFlag == 0)) printf("Using FFT to compute Yukawa interaction.\n");
  if ((myid == 0)&&(gFlag != 0)) printf("Using explicit formula for FT of Yukawa interaction.\n");

  if (myid != np-1) MPI_Send(&flag, 1, INT, myid+1, myid+1, COMM);
  MPI_Barrier(COMM);
  free(x); free(y); free(z);  
}
fftw_complex step(ctrl_ptr ctrl)
{
  if ((ctrl->heat) == 1) { 
    return (-1.I);
  } else return 1.;
}

